AIM : To Visualize the satellite products (SST, KPAR) as well as the in situ sites in Lyttelton harbour.

  1. KPAR is estimated using satellite measurements (MODIS AQUA) of absorption and backscattering coefficients. Their spatial resolution is 500m-pixel wide.
  2. SST data come from MODIS AQUA with a nominal resolution of 1km but has been downscalled to 500m-pixel wide.
  3. The in situ sites refer to cross-shore transects where ecological monitoring of presence/absence and abundance of macroalgae has been conducted by the Cawthron institute and related to the dredging of the harbour.

Document Note : The maps, plots and app within this document are interactive so make sure you give them a play like zooming in and out in the maps but also on the plots. Clicking on the legend allows to only select and display the time series needed.

Table of contents

  1. KPAR
  1. SST
  1. Scenario 1
  1. Scenario 2
  1. Scenario 1 vs Scenario 2

KPAR

The KPAR dataset consists of monthly means of the attenuation coefficient estimated using the absorption and backscattering coefficient using MODIS-AQUA measurements of the emergent flux, radiance leaving the water. The dataset runs from July 2002 to March 2019.

Mean

## Read site sheet
site <- read.csv(file=paste0(getwd(),'/Coords_sites.csv'))
##

## KPAR stack and mean
kpar_mean <- raster('A200207_201903_KPAR_Mean_Lyttelton_QAA.tif')
kpar <- stack('A200207_201903_KPAR_MO_Lyttelton_QAA.tif')
##

## Bathy of the zone NIWA grid
## Download Bathy, plot it with sites and transform as contour
bathy <- raster('Bathy_Lyttelton.tif')
bathy[bathy>0] <- NA
bathycontour <- rasterToContour(bathy)
bathycontour <- spTransform(bathycontour,crs("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0"))
##

## Plot KPAR mean and sites on leaflet map
pal <- colorNumeric(rev(rainbow(10)), values(kpar_mean),
                    na.color = "transparent")

m <- leaflet(site) %>% setView(lng = 173, lat = -43.55, zoom = 10) %>%
  addTiles()  %>%# Print the map
  addScaleBar(position = "bottomright",options = scaleBarOptions(imperial=F)) %>% 
  addMouseCoordinates() %>%
  addMarkers(site, lat = ~lat,lng = ~lon, popup = ~name) %>%
  addRasterImage(kpar_mean,layerId = "Kd", col=pal, opacity = 0.8,project=T) %>%
  addLegend(pal = pal, values = values(kpar_mean),title = "Kd (m-1)") %>%
  addImageQuery(x=kpar_mean,layerId = "Kd",type='click',project=T,position='bottomleft') %>%
  addPolylines(data=bathycontour,color = "black",opacity=0.2,popup = bathycontour$level)
m

Figure 1 -Mean KPAR product around Lyttelton harbour, NZ.

Key observations:

  • Some high values for pixels close to shore. This over-estimation could be due to the ‘land adjacency effect’ which is due to the land being an near infra-red (NIR) bright target which create a NIR halo affecting the atmospheric correction (Matt explanation).
  • To think about a method to take different pixes into account for the estimation of the KPAR value a the site. Maybe a compromise to find between the distance to the site and the availability of the pixel. Also possible to take more pixels into account.

Availability of pixels

kpar_NA_sum <- raster('A200207_201903_KPAR_PixAvailability_Lyttelton.tif')

kpar_NA_sum_proj <- projectRaster(kpar_NA_sum,crs=crs("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0"),method='ngb')

pal <- colorNumeric(rev(heat.colors(100)), values(kpar_NA_sum_proj),
                    na.color = "transparent")

m <- leaflet(site) %>% setView(lng = 173, lat = -43.55, zoom = 10) %>%
  addTiles()  %>%# Print the map
  addScaleBar(position = "bottomright",options = scaleBarOptions(imperial=F)) %>% 
  addMouseCoordinates() %>%
  addMarkers(site, lat = ~lat,lng = ~lon, popup = ~name) %>%
  addRasterImage(kpar_NA_sum_proj,layerId = "Number of pixels available", col=pal, opacity = 0.8) %>%
  addLegend(pal = pal, values = values(kpar_NA_sum_proj),title = "Number of pixels available") %>%
  addImageQuery(x=kpar_NA_sum_proj,layerId =  "Number of pixels available",type='click',position='bottomleft') %>%
  addPolylines(data=bathycontour,color = "black",opacity=0.2,popup = bathycontour$level)
m

Figure 2 - Availability of pixels for the KPAR product around Lyttelton harbour, NZ.

Key observations:

  • The “holes” offshore in the pixels are likely to be due because of the reprojection for leaflet. Nothing to worry about.
  • Fewer pixels available close to shore/sites. As said previous section, to think about a method that takes into account the availability of pixels and distance to site.

SST

The SST dataset runs from July 2002 to March 2019 and consists of monthly means. I don’t know from which satellite it is from but could retrieve that at the occasion.

Mean

## SST stack and mean 
sst <- stack('A200207_201903_SST_MO_Lyttelton.tif')
sst_mean <- raster('A200207_201903_SST_Mean_Lyttelton.tif')
##

## Plot SST mean and sites on leaflet map
pal <- colorNumeric(rev(rainbow(100)), values(sst_mean),
                    na.color = "transparent")

m <- leaflet(site) %>% setView(lng = 173, lat = -43.55, zoom = 10) %>%
  addTiles()  %>%# Print the map
  addScaleBar(position = "bottomright",options = scaleBarOptions(imperial=F)) %>% 
  addMouseCoordinates() %>%
  addMarkers(site, lat = ~lat,lng = ~lon, popup = ~name) %>%
  addRasterImage(sst_mean,layerId = "SST (degC)", col=pal, opacity = 0.8,project=T) %>%
  addLegend(pal = pal, values = values(sst_mean),title = "SST (degC)") %>%
  addImageQuery(x=sst_mean,layerId =  "SST (degC)",type='click',project=T,position='bottomleft') %>%
  addPolylines(data=bathycontour,color = "black",opacity=0.2,popup = bathycontour$level)
m

Figure 3 - Mean SST product around Lyttelton harbour, NZ.

Key observations:

  • Odd values in coastal pixels. Either too high? (in Lyttelton harbour and way out of bays) or too low? (within the bays).
  • As previously, to think about a method to get the SST estimated value by satellite to the sites.

Availability of pixels

sst_NA_sum <- raster('A200207_201903_SST_PixAvailability_Lyttelton.tif')

pal <- colorNumeric(rev(heat.colors(100)), values(sst_NA_sum),
                    na.color = "transparent")

m <- leaflet(site) %>% setView(lng = 173, lat = -43.55, zoom = 10) %>%
  addTiles()  %>%# Print the map
  addScaleBar(position = "bottomright",options = scaleBarOptions(imperial=F)) %>% 
  addMouseCoordinates() %>%
  addMarkers(site, lat = ~lat,lng = ~lon, popup = ~name) %>%
  addRasterImage(sst_NA_sum,layerId = "Number of pixels available", col=pal, opacity = 0.8,project=T) %>%
  addLegend(pal = pal, values = values(sst_NA_sum),title = "Number of pixels available") %>%
  addImageQuery(x=sst_NA_sum,layerId =  "Number of pixels available",type='click',position='bottomleft',project=T) %>%
  addPolylines(data=bathycontour,color = "black",opacity=0.2,popup = bathycontour$level)
m

Figure 4 - Availability of pixels for the SST product around Lyttelton harbour, NZ.

Key observations:

  • Less pixels available where dodgy values are suspected.

Conclusion and further Work

  • Because the sat data is available from 07/2002 to 03/2019, to define a time frame to refine the dataset. Use the surveys dates, maybe to focus on ecological data during the baseline period (19-20/01/17 BL1 to 5-8/12/17 BL3) and do another analysis with the data after dredging phase (DP1 only date to fell into sat data dataset range).
  • I could then regenerate monthly means within the new time frame and start exploring different scenarios for assigning a value to sites.

Different scenarios:

  • 1 : Extract values (KPAR and SST) at site pixel -> Will be not realistic because of missing values due to shore and effect of land.
  • 2 : Extract values at pixel wiht highest availability from cluster of 16 pixels adjacent to site location.
  • 2’ : Extract values at closest pixel offshore above threshold of certain pixel availability.
  • 3 : Extract values from cluster of pixels at certain distance of site and above threshold of pixel availability.

Scenario 1

In order to have estimates of KPAR and SST for each sites, we investigate here the scenario of extracting the time series of the pixel where the site is. This approach is supposed to fail as the sites are very coastal so prone to pixels unavailability due to cloud coverage as well as the land adjacency effect which leads to failure in the atmospheric correction process.

KPAR

## Work on Lat/Lon of site coordinates, reptroject to utm
# Reprojection of sites lon/lat into tmerc (EBED crs)
LatLong_site <- data.frame(Y=site$lat,X=site$lon)
names(LatLong_site) <- c("Y","X")
coordinates(LatLong_site) <- ~ X + Y # longitude first
proj4string(LatLong_site) <- crs("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0")
utm_sites <- spTransform(LatLong_site, crs( "+proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"))
##

## Extract values of pixel availability at sites location
kpar_NA_sum_velox <- velox('A200207_201903_KPAR_PixAvailability_Lyttelton.tif')
kpar_NA_sum_sites <- kpar_NA_sum_velox$extract_points(utm_sites)
sst_NA_sum_velox <- velox('A200207_201903_SST_PixAvailability_Lyttelton.tif')
sst_NA_sum_sites <- sst_NA_sum_velox$extract_points(utm_sites)
##

## KPAR
kpar_velox <- velox('A200207_201903_KPAR_MO_Lyttelton_QAA.tif')
kpar_sites <- kpar_velox$extract_points(utm_sites)
kpar_sites_df <- data.frame(t(kpar_sites))

colnames(kpar_sites_df) <- site$name
rownames(kpar_sites_df) <- NULL
dateseq <- seq.Date(as.Date("2002/7/1"), by = "month", length.out = dim(kpar_sites)[2])
kpar_sites_df$Date <- dateseq

kpar_sites_tb <- as_tibble(kpar_sites_df)
kpar_sites_tb <- kpar_sites_tb %>% pivot_longer(-Date,names_to="Site",values_to='Kpar') %>% group_by(Date)

q <- ggplot(kpar_sites_tb,aes(Date,Kpar,color=Site)) + geom_point() + geom_line() + labs(y="Kpar (/m)", x = "Date")
ggplotly(q)

Figure 5 - Time series of KPAR (/m) at the 21 sites around Lyttelton harbour, NZ. The values are extracted from the pixels where the sites are (Scenario 1).

Relevant observations:

  • No values for sites: PL14, PL03, PB10, PB11, PB03, PB02.
  • Very frequent gaps in time series.
  • Most of the values are over 0.2 /m and under 1 /m. Possible overestimation.
## KPAR mean, sd, median, min and max values over period
summary_kpar <- kpar_sites_tb %>% group_by(Site) %>% summarise(Mean=mean(Kpar,na.rm=T),Sd=sd(Kpar,na.rm=T),Median=median(Kpar,na.rm=T),Min=min(Kpar,na.rm=T),Max=max(Kpar,na.rm=T))
## `summarise()` ungrouping output (override with `.groups` argument)
summary_kpar$PixAv <- kpar_NA_sum_sites
kable(summary_kpar)
Site Mean Sd Median Min Max PixAv
BP01 0.4807346 0.1439982 0.4621158 0.2263836 0.9213209 86
BP02 0.5245741 0.1420317 0.5049014 0.2587270 0.9719090 83
BP03 0.5477195 0.1388396 0.5400223 0.3038899 0.9552073 48
BP05 0.3739728 0.1198098 0.3444114 0.2331560 0.7454379 47
BP06 0.4466732 0.1221743 0.4181812 0.2577299 0.9109438 106
BP08 0.4077619 0.1305797 0.3795969 0.1555093 1.0791025 90
BP10 0.4126234 0.1417000 0.3839348 0.2120700 0.8721717 67
BP13 0.4603269 0.1148118 0.4427282 0.2227825 0.7604557 90
BP14 0.6887431 0.2607180 0.6435227 0.2820656 1.4770980 70
LH01 0.5052427 0.1413214 0.5029341 0.2142351 0.8423380 73
LH02 0.4657920 0.1432052 0.4895723 0.2473924 0.6401853 9
LH07 0.6919909 0.2042824 0.6676622 0.3572536 1.4077549 60
LH10 0.6233862 0.1094229 0.6374094 0.4760661 0.7800066 7
PB02 NaN NA NA Inf -Inf 0
PB03 NaN NA NA Inf -Inf 0
PB10 NaN NA NA Inf -Inf 0
PB11 NaN NA NA Inf -Inf 0
PL02 0.8797934 0.2362799 0.9601306 0.6019558 1.1715964 5
PL03 NaN NA NA Inf -Inf 0
PL14 NaN NA NA Inf -Inf 0
PL16 0.8977092 0.2799484 0.9074726 0.5237424 1.2508827 8

SST

sst_velox <- velox('A200207_201903_SST_MO_Lyttelton.tif')
sst_sites <- sst_velox$extract_points(utm_sites)
sst_sites_df <- data.frame(t(sst_sites))

colnames(sst_sites_df) <- site$name
rownames(sst_sites_df) <- NULL
dateseq <- seq.Date(as.Date("2002/7/1"), by = "month", length.out = dim(sst_sites)[2])
sst_sites_df$Date <- dateseq

sst_sites_tb <- as_tibble(sst_sites_df)
sst_sites_tb <- sst_sites_tb %>% pivot_longer(-Date,names_to="Site",values_to='SST') %>% group_by(Date)

p <- ggplot(sst_sites_tb,aes(Date,SST,color=Site)) + geom_point() + geom_line() + labs(y="SST (degC)", x = "Date")
ggplotly(p)

Figure 6 - Time series of SST (degC) at the 21 sites around Lyttelton harbour, NZ. The values are extracted from the pixels where the sites are (Scenario 1).

Relevant observations:

  • No values for sites: PL03, PB11, PB10, PB03, PB02.
  • Very frequent gaps in time series. Some sites only have a few points (PL14), whereas others have only a few missing values (All BP for instance).
  • Most of the values are over 10 degC (Winter) and under 20 degC (Summer). Seems relatively good, however let’s have a look at the mean for every stations and see how the gaps in time series will affect the mean.
## SST mean, sd, median, min and max values over period
summary_sst <- sst_sites_tb %>% group_by(Site) %>% summarise(Mean=mean(SST,na.rm=T),Sd=sd(SST,na.rm=T),Median=median(SST,na.rm=T),Min=min(SST,na.rm=T),Max=max(SST,na.rm=T))
## `summarise()` ungrouping output (override with `.groups` argument)
summary_sst$PixAv <- sst_NA_sum_sites
kable(summary_sst)
Site Mean Sd Median Min Max PixAv
BP01 16.72645 3.505492 16.85250 8.695 24.035 100
BP02 16.93592 3.504716 17.20750 9.870 23.340 102
BP03 16.67018 3.449685 17.54500 8.170 23.110 51
BP05 14.50368 3.624410 14.19000 8.715 19.910 19
BP06 16.46508 3.629748 17.31250 8.520 23.965 118
BP08 15.74805 3.667408 16.24000 7.805 22.180 107
BP10 16.96604 3.810822 17.72750 8.300 23.660 78
BP13 14.59495 4.071016 14.70000 7.100 22.625 77
BP14 16.18449 4.400136 16.31474 9.260 25.740 50
LH01 14.94987 4.823432 15.33750 7.460 25.455 74
LH02 13.81568 5.620958 12.02000 5.640 24.265 26
LH07 18.14572 4.455413 18.60000 8.575 26.540 79
LH10 22.14909 3.711519 22.14909 12.875 30.245 23
PB02 NaN NA NA Inf -Inf 0
PB03 NaN NA NA Inf -Inf 0
PB10 NaN NA NA Inf -Inf 0
PB11 NaN NA NA Inf -Inf 0
PL02 17.39167 4.701594 16.73083 9.650 24.750 10
PL03 NaN NA NA Inf -Inf 0
PL14 7.54500 0.000000 7.54500 7.545 7.545 2
PL16 18.87106 4.213544 20.63000 9.630 24.130 27

Conclusion and further Work

  • The relative high number of missing values as well as the possible impact of the adjacent land lead to the conclusion that the scenario 1-way of getting estimated values for each site does not seem to show realistic values.
  • Let’s carry on the exploration and analyse the other scenarios.

Scenario 2

In order to have estimates of KPAR and SST for each sites, we investigate here the scenario of extracting the time series from a pixel with the maximum availability (using the pixel availability raster) within a cluster of 16 pixels adjacent to the site location (17 pixels including site, see what’s the structure of the cluster).

KPAR

## KPAR
ngb_pix_maxPixAv_index_save <- c()
sites_kpar_mean_save <- c()
for (i in 1:dim(site)[1]) {
  sites_tmerc <- data.frame(utm_sites)[i,]
  sites_kpar_mean <- raster::extract(kpar_mean,sites_tmerc,cellnumbers=T) #Get cell number of pixel where site
  ngb_pixels <- raster::adjacent(kpar_mean,cells=sites_kpar_mean[,1],pairs=F,directions=16,include=T)
  
  ngb_pix_maxPixAv_index <- ngb_pixels[which(kpar_NA_sum[ngb_pixels]==max(kpar_NA_sum[ngb_pixels]))] #Return index of pixel in cluster of 16 pixels around that have max pixel availability
  # Issues with PB10 and PB3 -> Seems like no pixels available in vicinity
  if (length(ngb_pix_maxPixAv_index)>1) { 
    #Manually solve the issue affecting original index to pixel
    ngb_pix_maxPixAv_index <- sites_kpar_mean[,1]
  }
  ngb_pix_maxPixAv_index_save <- c(ngb_pix_maxPixAv_index_save,ngb_pix_maxPixAv_index)
}

# Position of the "new" pixels where extractions of scenario 2 is from
pix_index_maxPixAv <- xyFromCell(kpar_mean,cell=ngb_pix_maxPixAv_index_save,spatial=F)#Convert to sp feature, to use velox::
Tmerc_site_s2 <- data.frame(Y=pix_index_maxPixAv[,2],X=pix_index_maxPixAv[,1])
coordinates(Tmerc_site_s2) <- ~ X + Y # longitude first
proj4string(Tmerc_site_s2) <- crs("+proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
LatLong_sites_s2 <- spTransform(Tmerc_site_s2,crs("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0"))
LatLong_sites_s2$name <- paste0(site$name,"_s2")
LatLong_sites_s2_df <- data.frame(LatLong_sites_s2)

# Extract time series from "new" pixels
pix_index_maxPixAv_sp <- xyFromCell(kpar_mean,cell=ngb_pix_maxPixAv_index_save,spatial=T)#Convert to sp feature, to use velox::
kpar_sites_maxPixAv <- kpar_velox$extract_points(pix_index_maxPixAv_sp)#Classic extraction process
kpar_sites_maxPixAv_df <- data.frame(t(kpar_sites_maxPixAv))

colnames(kpar_sites_maxPixAv_df) <- LatLong_sites_s2$name#paste0(site$name[1],'_maxPixAv')
rownames(kpar_sites_maxPixAv_df) <- NULL
dateseq <- seq.Date(as.Date("2002/7/1"), by = "month", length.out = dim(kpar_sites_maxPixAv_df)[1])
kpar_sites_maxPixAv_df$Date <- dateseq

kpar_sites_maxPixAv_tb <- as_tibble(kpar_sites_maxPixAv_df)
kpar_sites_maxPixAv_tb <- kpar_sites_maxPixAv_tb %>% pivot_longer(-Date,names_to="Site",values_to='Kpar') %>% group_by(Date)

q <- ggplot(kpar_sites_maxPixAv_tb,aes(Date,Kpar,color=Site)) + geom_point() + geom_line() + labs(y="Kpar (/m)", x = "Date")
ggplotly(q)

Figure 7 - Time series of KPAR (/m) at the 21 sites around Lyttelton harbour, NZ. The values are extracted from the pixels with the highest availability in a cluster of 16 pixels adjacent to each sites. Scenario 2.

m <- leaflet(data=LatLong_sites_s2_df) %>% setView(lng = 173, lat = -43.55, zoom = 10) %>%
  addTiles()  %>%# Print the map
  addMarkers(LatLong_sites_s2_df, lat = ~Y,lng = ~X, popup = ~name)
m

Figure 8 - Position of the Pixels where time series are extracted in the scenario 2 case.

Relevant observations:

  • TO FILL (Less gap, Position of the new pixels, PB03 and PB10 no value coz no ngb)
## KPAR mean, sd, median, min and max values over period
summary_kpar_s2 <- kpar_sites_maxPixAv_tb %>% group_by(Site) %>% summarise(Mean=mean(Kpar,na.rm=T),Sd=sd(Kpar,na.rm=T),Median=median(Kpar,na.rm=T),Min=min(Kpar,na.rm=T),Max=max(Kpar,na.rm=T))
## `summarise()` ungrouping output (override with `.groups` argument)
summary_kpar_s2$PixAv <- kpar_NA_sum[ngb_pix_maxPixAv_index_save]
kable(summary_kpar_s2)
Site Mean Sd Median Min Max PixAv
BP01_s2 0.3736091 0.0757962 0.3536769 0.2344388 0.7235438 168
BP02_s2 0.3566422 0.0690696 0.3449551 0.2511407 0.7337249 184
BP03_s2 0.3382366 0.0698441 0.3276764 0.2026570 0.6407034 169
BP05_s2 0.3479043 0.0756010 0.3367585 0.2061284 0.6099731 163
BP06_s2 0.3381993 0.0727750 0.3218081 0.2134089 0.5809690 190
BP08_s2 0.3111466 0.0691485 0.2943282 0.1995939 0.6331935 190
BP10_s2 0.3616482 0.0944326 0.3516572 0.1610703 0.9613709 185
BP13_s2 0.3749535 0.0785294 0.3640774 0.1974660 0.6571044 180
BP14_s2 0.3566313 0.0928171 0.3355087 0.1544371 0.6814184 177
LH01_s2 0.4792453 0.1089673 0.4638782 0.2742959 0.8294032 138
LH02_s2 0.5040968 0.1202039 0.4734725 0.2543523 0.9888393 145
LH07_s2 0.5214943 0.1296387 0.4985805 0.2543523 0.9773747 146
LH10_s2 0.4840434 0.1156862 0.4757362 0.2742959 0.8221809 136
PB02_s2 0.4004347 0.0580121 0.4073745 0.3243970 0.4625929 4
PB03_s2 NaN NA NA Inf -Inf 0
PB10_s2 NaN NA NA Inf -Inf 0
PB11_s2 1.1407515 NA 1.1407515 1.1407515 1.1407515 1
PL02_s2 0.6564987 0.2188256 0.6817867 0.2699000 1.0174295 41
PL03_s2 1.0310847 0.2605226 1.1485419 0.6019558 1.3445294 11
PL14_s2 1.0310847 0.2605226 1.1485419 0.6019558 1.3445294 11
PL16_s2 0.4569534 0.1356378 0.4164062 0.2840720 0.9177421 87

SST

## SST
ngb_pix_maxPixAv_index_save <- c()
sites_sst_mean_save <- c()
for (i in 1:dim(site)[1]) {
  sites_tmerc <- data.frame(utm_sites)[i,]
  sites_sst_mean <- raster::extract(kpar_mean,sites_tmerc,cellnumbers=T) #Get cell number of pixel where site
  ngb_pixels <- raster::adjacent(sst_mean,cells=sites_sst_mean[,1],pairs=F,directions=16,include=T)
  
  ngb_pix_maxPixAv_index <- ngb_pixels[which(sst_NA_sum[ngb_pixels]==max(sst_NA_sum[ngb_pixels]))] #Return index of pixel in cluster of 16 pixels around that have max pixel availability
  # Issues with PB10 and PB3 -> Seems like no pixels available in vicinity
  if (length(ngb_pix_maxPixAv_index)>1) { 
    #Manually solve the issue affecting original index to pixel
    ngb_pix_maxPixAv_index <- sites_sst_mean[,1]
  }
  ngb_pix_maxPixAv_index_save <- c(ngb_pix_maxPixAv_index_save,ngb_pix_maxPixAv_index)
}

# Position of the "new" pixels where extractions of scenario 2 is from
pix_index_maxPixAv_sst <- xyFromCell(sst_mean,cell=ngb_pix_maxPixAv_index_save,spatial=F)#Convert to sp feature, to use velox::
Tmerc_site_s2_sst <- data.frame(Y=pix_index_maxPixAv_sst[,2],X=pix_index_maxPixAv_sst[,1])
coordinates(Tmerc_site_s2_sst) <- ~ X + Y # longitude first
proj4string(Tmerc_site_s2_sst) <- crs("+proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
LatLong_sites_s2_sst <- spTransform(Tmerc_site_s2_sst,crs("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0"))
LatLong_sites_s2_sst$name <- paste0(site$name,"_s2")
LatLong_sites_s2_sst_df <- data.frame(LatLong_sites_s2_sst)

# Extract time series from "new" pixels
pix_index_maxPixAv_sp <- xyFromCell(sst_mean,cell=ngb_pix_maxPixAv_index_save,spatial=T)#Convert to sp feature, to use velox::
sst_sites_maxPixAv <- sst_velox$extract_points(pix_index_maxPixAv_sp)#Classic extraction process
sst_sites_maxPixAv_df <- data.frame(t(sst_sites_maxPixAv))

colnames(sst_sites_maxPixAv_df) <- LatLong_sites_s2_sst$name#paste0(site$name[1],'_maxPixAv')
rownames(sst_sites_maxPixAv_df) <- NULL
dateseq <- seq.Date(as.Date("2002/7/1"), by = "month", length.out = dim(sst_sites_maxPixAv_df)[1])
sst_sites_maxPixAv_df$Date <- dateseq

sst_sites_maxPixAv_tb <- as_tibble(sst_sites_maxPixAv_df)
sst_sites_maxPixAv_tb <- sst_sites_maxPixAv_tb %>% pivot_longer(-Date,names_to="Site",values_to='SST') %>% group_by(Date)

q <- ggplot(sst_sites_maxPixAv_tb,aes(Date,SST,color=Site)) + geom_point() + geom_line() + labs(y="SST (degC)", x = "Date")
ggplotly(q)

Figure 9 - Time series of SST (degC) at the 21 sites around Lyttelton harbour, NZ. The values are extracted from the pixels with the highest availability in a cluster of 16 pixels adjacent to each sites. Scenario 2.

m <- leaflet(data=LatLong_sites_s2_sst_df) %>% setView(lng = 173, lat = -43.55, zoom = 10) %>%
  addTiles()  %>%
  addMarkers(LatLong_sites_s2_sst_df, lat = ~Y,lng = ~X, popup = ~name)
m

Figure 10 - Position of the pixels where time series of SST are extracted from in the scenario 2 case.

Relevant observations:

  • TO FILL (Less gap, Position of the new pixels, PB03 and PB10 no value coz no ngb, PL02 new pixel is in LH.)
## KPAR mean, sd, median, min and max values over period
# SST mean, sd, median, min and max values over period
summary_sst_s2 <- sst_sites_maxPixAv_tb %>% group_by(Site) %>% summarise(Mean=mean(SST,na.rm=T),Sd=sd(SST,na.rm=T),Median=median(SST,na.rm=T),Min=min(SST,na.rm=T),Max=max(SST,na.rm=T))
## `summarise()` ungrouping output (override with `.groups` argument)
summary_sst_s2$PixAv <- sst_NA_sum[ngb_pix_maxPixAv_index_save]
kable(summary_sst_s2)
Site Mean Sd Median Min Max PixAv
BP01_s2 14.55268 3.2726353 15.1200 7.785000 20.520 172
BP02_s2 16.93592 3.5047162 17.2075 9.870000 23.340 102
BP03_s2 14.35078 3.2165821 14.6975 8.120000 20.775 166
BP05_s2 14.69212 3.3178140 14.9900 8.479999 21.245 142
BP06_s2 13.86971 3.3077064 13.9600 8.264999 19.750 194
BP08_s2 13.94149 3.2028660 14.2700 8.160000 19.295 187
BP10_s2 13.93659 3.1511820 14.1350 8.255000 19.590 184
BP13_s2 14.47864 3.5653955 14.6300 8.264999 21.900 160
BP14_s2 14.44132 3.4180549 14.6200 8.040000 21.285 157
LH01_s2 14.78414 3.8983909 15.4800 7.535000 22.210 153
LH02_s2 14.93382 4.3012935 15.1625 7.610000 25.065 164
LH07_s2 14.82535 3.7880783 15.3200 7.580000 22.195 155
LH10_s2 14.95026 3.8885392 15.6750 7.770000 22.660 151
PB02_s2 13.69000 0.3400002 13.6900 13.349999 14.030 3
PB03_s2 NaN NA NA Inf -Inf 0
PB10_s2 NaN NA NA Inf -Inf 0
PB11_s2 21.83500 6.1904860 22.7275 13.520000 28.365 4
PL02_s2 19.21261 4.1895382 19.9050 9.525000 28.535 58
PL03_s2 19.00000 3.9163583 18.7625 13.145000 25.945 14
PL14_s2 17.37219 5.2478802 17.8000 5.345000 25.330 17
PL16_s2 16.41398 3.6061747 16.9075 8.139999 22.940 106

Scenario 1 vs Scenario 2

The following tables show the difference of the values from the scenario 1 minus scenario 2 for each sites. Negative values mean that the scenario 2 leads to an increase in the parameter/value from scenario 1.

KPAR

## KPAR - Relative difference: s1 - s2
is.na(summary_kpar) <- sapply(summary_kpar, is.infinite)#Replace INF value by NA
is.na(summary_kpar_s2) <- sapply(summary_kpar_s2, is.infinite)#Replace INF value by NA
rdiff_s1s2_kpar <- summary_kpar[,2:7] - summary_kpar_s2[,2:7]
rdiff_s1s2_kpar$Site <- summary_kpar$Site
kable(rdiff_s1s2_kpar)
Mean Sd Median Min Max PixAv Site
0.1071255 0.0682020 0.1084389 -0.0080553 0.1977770 -82 BP01
0.1679319 0.0729621 0.1599463 0.0075863 0.2381842 -101 BP02
0.2094828 0.0689955 0.2123459 0.1012329 0.3145038 -121 BP03
0.0260685 0.0442089 0.0076529 0.0270276 0.1354648 -116 BP05
0.1084739 0.0493994 0.0963730 0.0443210 0.3299748 -84 BP06
0.0966152 0.0614312 0.0852687 -0.0440846 0.4459090 -100 BP08
0.0509752 0.0472675 0.0322776 0.0509997 -0.0891992 -118 BP10
0.0853734 0.0362824 0.0786507 0.0253166 0.1033512 -90 BP13
0.3321118 0.1679009 0.3080139 0.1276285 0.7956796 -107 BP14
0.0259974 0.0323541 0.0390559 -0.0600609 0.0129349 -65 LH01
-0.0383048 0.0230013 0.0160998 -0.0069599 -0.3486540 -136 LH02
0.1704966 0.0746438 0.1690816 0.1029013 0.4303802 -86 LH07
0.1393428 -0.0062633 0.1616731 0.2017702 -0.0421743 -129 LH10
NaN NA NA NA NA -4 PB02
NaN NA NA NA NA 0 PB03
NaN NA NA NA NA 0 PB10
NaN NA NA NA NA -1 PB11
0.2232947 0.0174544 0.2783440 0.3320558 0.1541669 -36 PL02
NaN NA NA NA NA -11 PL03
NaN NA NA NA NA -11 PL14
0.4407557 0.1443106 0.4910663 0.2396705 0.3331407 -79 PL16
##

Key observations:

  • Considerable gain of pixel availability, except for 6 sites.
  • Not relevant methods for 6 sites where originally NA.
  • Scenario 2 leads to decrease of KPAR mean.

SST

## SST - Relative difference: s1 - s2
is.na(summary_sst) <- sapply(summary_sst, is.infinite)#Replace INF value by NA
is.na(summary_sst_s2) <- sapply(summary_sst_s2, is.infinite)#Replace INF value by NA
rdiff_s1s2_sst <- summary_sst[,2:7] - summary_sst_s2[,2:7]
rdiff_s1s2_sst$Site <- summary_sst$Site
kable(rdiff_s1s2_sst)
Mean Sd Median Min Max PixAv Site
2.1737668 0.2328568 1.7325006 0.9099998 3.5149994 -72 BP01
0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0 BP02
2.3194007 0.2331030 2.8475003 0.0500002 2.3349991 -115 BP03
-0.1884315 0.3065963 -0.8000002 0.2350006 -1.3349991 -123 BP05
2.5953728 0.3220413 3.3524995 0.2550001 4.2150002 -76 BP06
1.8065653 0.4645421 1.9700003 -0.3550000 2.8850002 -80 BP08
3.0294443 0.6596399 3.5924993 0.0450001 4.0699997 -106 BP10
0.1163033 0.5056204 0.0700002 -1.1649995 0.7250004 -83 BP13
1.7431729 0.9820810 1.6947432 1.2200003 4.4549999 -107 BP14
0.1657225 0.9250412 -0.1424999 -0.0749998 3.2450008 -79 LH01
-1.1181394 1.3196644 -3.1424999 -1.9699998 -0.7999992 -138 LH02
3.3203697 0.6673344 3.2800007 0.9949999 4.3449993 -76 LH07
7.1988307 -0.1770199 6.4740906 5.1050000 7.5849991 -128 LH10
NaN NA NA NA NA -3 PB02
NaN NA NA NA NA 0 PB03
NaN NA NA NA NA 0 PB10
NaN NA NA NA NA -4 PB11
-1.8209427 0.5120556 -3.1741657 0.1250000 -3.7849998 -48 PL02
NaN NA NA NA NA -14 PL03
-9.8271873 -5.2478802 -10.2549996 2.1999998 -17.7850003 -15 PL14
2.4570789 0.6073695 3.7224998 1.4900007 1.1900005 -79 PL16

Key observations:

  • Considerable gain of pixel availability, except for 5 sites.
  • Not relevant methods for 5 sites where originally NA
  • Scenario 2 leads to decrease of SST mean except at 4 sites (BP05, LH02, PL02, PL14).
  • PB02 wtf, to check